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ABSTRACT 

This paper presents a novel semi-automatic image processing technique to estimate accu- 
rately, and objectively, the disc parameters of a planetary body on an astronomical image. The 
method relies on the detection of the limb and/or the terminator of the planetary body with 
the VOronoi Image SEgmentation (VOISE) algorithm ( |Guio, P. and Achilleos , N. 2009). The 
resulting map of the segmentation is then used to identify the visible boundary of the plane- 
tary disc. The segments comprising this boundary are then used to perform a "best" fit to an 
algebraic expression for the limb and/or terminator of the body. We find that we are able to 
locate the centre of the planetary disc with an accuracy of a few tens of one pixel. The method 
thus represents a useful processing stage for auroral "imaging" based studies. 

Key words: methods: data analysis — methods: miscellaneous — methods: statistical — 
techniques: image processing. 



1 INTRODUCTION 

During the last two decades, the Hubble Space Telescope (HST) 
has provided resolved images of both Jupiter and Saturn in the ul- 
traviolet (UV) spectral region. Such images capture with high sen- 
sitivity and high resolution, the spectacular auroral phenomena oc- 
curring in the polar regions of the gas giants as a result of energetic 
magnetospheric particles raining down onto the planet's upper at- 
mosphere. Auroral images have become a particularly useful diag- 
nostic tool for morphological characterisations of the aurora and 
its boundaries. This is a crucial prerequisite for identifying the au- 
rora's physical origin (e.g. Prange et al. 1996, 1998; Grodent et al. 
|2003b|a[|CTarke et al.|2005||Badman et al.|2008||Lamy et al.|200'9l 
Imaging is also complementary to in situ measurements of the 
plasma environment provided, e.g. by the Cassini spacecraft, cur- 
rently orbiting Saturn. Combining remote imaging with in situ data 
allows the study of magnetospheric processes and how they affect 
the planet's upper atmosphere, and ionosphere via the planet's mag- 
netic field (Dougherty et al.|1998||Clarke et al.|2002"l|Bunce et al.| 
|2008| |Talboys et al.||2009), a nd the footprint auroral emission of 
satellites (Clarke et al.|2002|[Bon"fond et al.|2007l|Wannawichian| 
|et al.||20"08| >. Such studies require accurate projection of the geo- 
graphic and geomagnetic coordinate systems of the planet onto the 
plane of the two-dimensional image. Auroral dynamics can be stud- 
ied using time series of images. For these purposes, the location of 
the planet centre needs to be known accurately. Unfortunately, HST 
pointing parameters are not generally known with sufficient accu- 
racy for such applications. The precision of the guide star catalogue 
together with the uncertainty in the start time of the tracking motion 
is on the order of 1 arc sec while it is desired to have an accuracy of 
the order of 1 pixel, i.e. 0.02-0.03 arc sec for the Space Telescope 

* E-mail:p. guio@ucl.ac.uk 



imaging spectrograph (STIS) and Advanced Camera for Surveys 
(ACS) instruments, in order to to locate any structure accurately or 
to build polar projections of the auroral emissions. 

In addition, ground-based observations with telescopes such 
as the NASA Infrared Telescope Facility (IRTF) and United King- 
dom Infra-Red Telescope (UKIRT), both located at Mauna Kea, 
Hawaii, have provided images of Jupiter and Saturn with resolved 
auroral structures in the infrared (IR) waveband. IR images and 
spectra also allow the study of the dynamics and morphology of 
the H^ molecular ion, a principal component of giant ionospheres 
(Miller et al. 2006). Again, the location of the planet centre needs 
to be known accurately to make use of these images, but for similar 
reasons as the HST case, the pointing parameters are not known 
with sufficient accuracy for the images from IRTF and UKIRT 
telescopes. The resolution of the IRTF imaging facility and the 
UKIRT imager-spectrometer (UIST) are respectively of the order 
of 0.04 arc sec pixel" 1 and 0.12 arc sec pixel - or better. 

The problem of the location of the planet on auroral images 
has been addressed by various authors and studies (e.g. IBonfond] 
|et~aLl|2007l |Nichols et aLl[2008] |Bonfond et aTpOOl?) but to our 
knowledge no published work provides any detailed description of 
the method used. Here we propose a novel semi-automatic method 
to estimate accurately and objectively the position, size and orien- 
tation of a planetary body. The method consists of three phases: 
(i) detection of the limb of the planet disc using our image seg- 
mentation algorithm VOronoi Image SEgmentation (VOISE) (see 
|Guio, R and Achi lleos, N. ( 2009 1 for details), (ii) selection of points 
(Voronoi seeds) from the VOISE map that surround the limb, and 
(iii) nonlinear fit (Levenberg-Marquardt algorithm) of the selected 
set of data from VOISE to a disc model. Phase (i) is performed once 
while phases (ii) and (iii) can be repeated in order to improve the 
accuracy. 

In section [2] we give analytic expressions for the projection 
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Figure 1. Sketch of the geometry of the planet and the observer. The eccen- 
tricity of the planetary ellipsoid is exaggerated for clarity. 



of the limb and the terminator in the sky-plane. In section [3] the 
method is developed. In section [4] we illustrate our method on IR 
images of Jupiter collected with the IRTF and UKIRT telescopes. 
We discuss the performance of our method and summarise our con- 
clusions in section[5] 



2 LIMB AND TERMINATOR EQUATIONS 

A planet's pressure surface can be modelled by an ellipsoid, more 
precisely an oblate spheroid, with equatorial radius (semi-major 
axis) r e and polar radius (semi-minor axis) r v , where r 2 , = 
r 2 (l— e 2 ), and e is the eccentricity of the spheroid. The parameters 
are readily available, for instance, from the NASA Navigation and 
Ancillary Information Facility SPICE system ( Acton 1996}. The 
planet rotation vector is assumed to be along the z-axis with posi- 
tive angular velocity oj as shown in Fig.[T] Without loss of general- 
ity, we can further assume that the observer is located in the (x, z)- 
plane, i.e. setting the longitude of the observer A©=0, and the ob- 
serving direction as the vector <5® = (cos 0, sin /}$) where /3® 
is the planetocentric latitude (sub-Earth latitude), the (negative) an- 
gle between the i-axis and J® in Fig. [T] The limb of the planet 
consists of the points on the planet surface where the local normal 
ns is perpendicular to the observing direction 5® , i.e. ns -5® =0. 
The entire limb is contained in a single plane and is an ellipse, and 
the normal vector to the plane containing the limb n l has coordi- 
nates (cos fi n , 0, sin P n ) where 



(l-e 2 )cos/3 a 



^(l_ e 2)2 cos 2 p^+sm 2 fa 



sin/3^ 



sin /3 a 



V(l-e 2 ) 2 cos 2 /3 +sm 2 /3 a 
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The ellipse of the limb can then be projected onto the plane of the 
sky [x 3 ,y s ), i-e. on a plane perpendicular to the observing direction 
(5© (see Fig.[T}. The projection of the limb in the sky-plane is also 
an ellipse with the following algebraic equation in the sky plane 
coordinate system 



where the j/ s -axis is chosen to lie in the (a;, z)-plane. It can be seen 
from Eq. |3j that the limb always appears with a semi-major axis 
equal to the equatorial radius of the planet and a semi-minor axis 
between the polar radius of the planet (in the case of an equatorial 
view f3(D — 0) and the equatorial radius of the planet (in the case 
of a polar view /3® = ±7r/2). Equivalently the eccentricity of the 
ellipse formed by the limb in the sky is et = e cos /3®. 

The terminator (boundary between day and night side) can be 
visualised as the limb for the direction corresponding to the loca- 
tion of the Sun Sq (i.e. such that the local normal fir is perpen- 
dicular to Sq) but projected onto the sky plane along the observ- 
ing direction <5®, i.e. directed towards Earth. The Sun direction is 
defined by its planetocentric latitude /3 (sub-solar latitude) and 
its relative longitude (solar phase angle) to the observing direction 
AA = Aq — A®. In the coordinate system (x',y',z') where z' is 
also aligned to the rotation axis of the planet but the plane (x' , z') 
has been rotated about the planet rotation axis to contain the Sun 
direction, the vector fix has coordinates (cos/3^ T , 0, sin/3^^) 
where 



cos j3. 



(1— e 2 ) cos/30 



nT v /(l-e 2 ) 2 cos 2 /3 +sin 2 /3 ' 
sin^ 



v /(l-e 2 ) 2 cos 2 /3 +sin 2 /3 



(4) 
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In the situation where 5® x fir = 0, the limb and the terminator 
are coincident. Otherwise the vector <5 x Aj is contained in the 
sky plane (since it is perpendicular to 5®) and in the plane of the 
terminator. Therefore the projections onto the sky-plane of the limb 
and the terminator intersect at two points called the cusps, and the 
line joining the two cusps has direction <5© x fir- 

The two cusp points define the major axis of the ellipse formed 
by the sky projection of the terminator. This shape is a tilted ellipse 
with tilt angle 9t with respect to the a; s -axis given by 



tan 



(1—e 2 ) cos Pq sin AA 
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The semi-major and semi-minor axes of the projection of the full 
terminator (i.e. its visible and invisible parts) onto the sky-plane are 



2 / 2 i 2\ .2 
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br = {u 2 +v 2 )tl. 



(7) 
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where the vector u — (u, v) lies in the direction defined by 

d xn T , 

u = (1—e 2 ) cos /3q cos AA sin /3® — sin /3 cos /3®, (9) 
v = (l-e 2 )cos/? sinAA, (10) 
and ti and ti are scalars analytically derivable from u and v 

r 2 (l-e 2 cos 2 /? ffi ) 



ti 



m 2 (1— e 2 cos 2 /3 ) + v 2 ' 

r 2 (l-el e )(ad-bc) 2 
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(11) 
(12) 



where e^ Q is the eccentricity of the ellipse formed by the termi- 
nator in its own plane and can be expressed as function of the Sun 
planetocentric latitude /3 and the eccentricity of the spheroid e 
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and where 

a = cosAA, (14) 

b = -sin/3ft T sinAA, (15) 

c = sin AAsin/3®, (16) 

d = sin/3^ r cos AA sin /3® + cos Pfy T cos/3 e (17) 

Finally the signed distance (as measured along fir) between 
a point (xt, Ht) of the projection of the terminator onto the sky- 
plane, and the plane of the terminator itself is given by 

Dt = (cos/?j^ T sin AX)xt+ 

(-cos/3^ t cos AAsin/3®+sin/3 AT cos/3 e )j/T, (18) 

and points with Dt>0 belong to the visible terminator (from the 
Earth observer's point of view) while points with Dt<0 are hid- 
den, and the case Dt~0 corresponds to cusp points. Similarly the 
signed distance (measured along ^0) between a point of the limb's 
projection onto the sky-plane (xl,Vl), and the plane perpendicular 
to the direction of the Sun Sq is given by 

Dl = (cos Sq sinAA)a;_L + 

(- cos 6 e cos A A sin /3 ffi + sin S Q C °^^f ^ Dl, (19) 

and points such that Dl>0 belong to the illuminated limb while 
points such that Dl<0 are in the shade, and the case Dl~0 corre- 
sponds to cusp points. 



3 PLANETARY DISC EXTRACTION METHOD 

As pointed out the proposed method to extract the orientation and 
shape of the planetary disc from an image consists of three phases, 
described in the following sections. 



3.1 Phase (i) VOISE image reduction 

The first stage consists of partitioning the image into regions, i.e. 
simplify and/or change the representation of an image into some- 
thing that is more physically meaningful and easier to analyse. 
VOISE is a dynamic algorithm for partitioning the underlying pixel 
grid of an image into regions according to a prescribed homogene- 
ity criterion ( |Guio, P. a nd Achilleos, N. 2009 1. A VOISE segmen- 
tation returns a map of the image in the form of a Voronoi dia- 
gram (VD) where each Voronoi region (VR) is a polygon, within 
which the data are homogeneous with respect to prescribed crite- 
ria. When running the VOISE segmentation algorithm on an image 
of a planetary object we expect that the transition region between 
the illuminated planet and the sky (i.e. the limb or the terminator) 
consists of a ring of relatively small Voronoi polygons, indicating 
that at this region, the intensity is changing very quickly over small 
spatial scales. In this representation we can classify the "ring" of 
tiny Voronoi polygons surrounding the larger central polygons as 
a cluster in itself that can be used for fitting a terminator and/or a 
limb. The map generated at the end of the VOISE division phase 
provides the largest number of seeds and smallest Voronoi poly- 
gons (see |Guio, P. and Achilleos7N^p009) for more detail). 



3.2 Phase (ii) points selection 

The selection of seeds from the computed VOISE map requires 
"crude" estimates for the planet centre (x c ,y c ), the equatorial ra- 
dius (semi-major axis) r e , the polar radius (semi-minor axis) r p 
and the tilt angle a. The seeds from the segmentation are consid- 
ered part of the neighbourhood of the limb and/or terminator if they 
lie inside a prescribed elliptic torus. A point belongs to the torus if 
its coordinates (x, y) fulfil the following inequalities 
/2 ,2 

a 2 b 2 

where e m and em (with e m < 1 < em) represent the inner and 
outer ellipses of the torus, and where (x',y') are obtained from 
(x, y) by translation with —(x c , y c ) followed by rotation with —a, 
i.e. 
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As well as satisfying the condition given by Eqs. l |20f|21[ l, the poly- 
gons used in the fitting procedure must also have a surface area 
smaller than a prescribed value (equivalently a "length scale" C 
smaller than a maximum prescribed value Cm)- The maximum 
length scale Cm has to be larger than the minimum distance be- 
tween seeds d m that is set for the VOISE division phase. 

It is also possible to filter out "bands" of polar angle in order 
to remove those regions with auroral emission clearly outside the 
planetary limb as illustrated in section]?] 

3.3 Phase (iii) fitting of points 

Fitting of quadrics (such as circles and ellipses) to a given set of 
points in the plane is a problem that arises in many application ar- 
eas, e.g. computer graphics, pattern recognition, coordinate meteo- 
rology. Many algorithms minimise a quantity in some least-square 
sense. Such fitting algorithms for quadrics can be separated into 
categories of "best fit" ("geometric fit") and "algebraic fit" ( |Gan-| 
|der et al.|1994| |Fitzgibbon et al.|1999") . In addition the clustering 
technique is another technique to fit an ellipse, such as methods 
based on the Hough transform (Yuen et aT[ l989 1. 

In the "best fit", the quantity to minimise is the geometric dis- 
tance between the fitted curve and the given set of points. In this 
case, curves may be represented in parametric form, which is well 
suited for minimising the sum of the squares of the distances. 

In the "algebraic fit" the curve is represented algebraically, 
i.e. in the plane by an equation of the form F(x, y) — 0. If a 
point is on the curve, then its coordinates (x, y) are a zero of the 
function F and represent an algebraic distance. These methods are 
usually equivalent to solving a linear system of equations subject 
to some constraint on the quadric coefficients (Bookstein 1979; 
|Taubin 1991 ; Fitzgib bon et al.| 19 99 1. The constraint may be such 
that the optimal solution is computed directly, and no iterations are 
needed. The disadvantage of the "algebraic fit" is that we are un- 
certain what we are minimising in a geometrical sense and in many 
cases those constraints lead to fits which are not invariant under Eu- 
clidean transformations such as translations and rotations, i.e. dif- 
ferent coordinate systems produce different fitted curves. Another 
limitation is that these methods can fit only one primitive (or shape) 
at a time, therefore the data should be segmented into a set of ba- 
sic shapes before fitting each shape independently. Nonetheless the 
algebraic solution is useful as an initial guess for the geometrical 
fit. 

We have developed a "best fit" tool based on the the 
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Levenberg-Marquardt method (Marquardt 1963), for solving non- 
linear least-square problems. This method performs a minimisa- 
tion of the sum of the squares of the weighted distances between 
the m selected seeds Si (i=l . . . m) from the Voronoi map, and 
the "best"-fitting curve with parametric representation (x, y) — 
f(</>;p) ( Bard|1974 Gander et al.||1994 . The minimisation con- 
sists in adjusting iteratively the set of curve parameters {4>i} (that 
locate the "best" points on the curve), together with the vector of 
global parameters p (that describe the global shape of the curve). 
Mathematically the function to minimise is written 



Q(4>i,4>2 



f(4HiP)\f 



(22) 



The parametric representation / for a circle and an ellipse are given 
by Eq. \23\ and Eq. \24\ respectively, of represents the variance 
of the location of a given seed. An estimate for this uncertainty 
can be readily computed as the mean distance from the seed to all 
the points within the VR. Alternatively a length scale d of the 
polygon can be inferred from the square root of its surface area 
( |Guio, P. and Achilleos, N.|20 09 1, and can be thought of as the "av- 
erage" section length through the polygon in all directions. Thus 
considering the disc with same surface area A as the polygon, the 
variance in distance of the disc points from its centre is given by 
a 2 = A/(2tv) = C 2 /2 (where surface area is used as the weight- 
ing factor for variance). This expression can thus be used as a rea- 
sonable approximation for the variance of the location of the seeds. 

The Levenberg-Marquardt method is optimised to switch con- 
tinuously from a method which quickly approaches the minimum 
(the steepest descent method), when far from the minimum, to a 
more precise but slower method (the Newton method), when ap- 
proaching the minimum. 

Finally we note that a priori knowledge about any of the pa- 
rameters can be used to constrain the fitting, otherwise an appropri- 
ate number of free parameters may be simultaneously determined 
from the fitting procedure. 

3.3.1 Fitting a circle 

The parametric form used for the circle is given by 
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where x c , y c and r are respectively the coordinates of the centre and 
the radius of the circle. Note that the values of <j> for each data point 
Si are updated along with the global parameters p = [x c ,y c , r] 
at each iteration. The parameter (j> represents the polar angle, mea- 
sured from the x-axis, of the line joining the centre (x c , y c ) of the 
circle to the point on the circle which is associated with the relevant 
data point. 



3.3.2 Fitting an ellipse 

The parametric form used for the ellipse is given by 

f{4>\ [x c ,y c ,a,b,a]) = 
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where x c , y c , a, b and a are respectively the coordinates of the 
centre, the semi-major axis, semi-minor axis and the tilt angle of 
the ellipse (angle measured from the s-axis to the semi-major axis). 
In this case the parameter </> does not represent the polar angle, 



measured from the x-axis, of the line joining the centre (x c , y c ) of 
the ellipse to the point on the ellipse. The parameter <fi is sometimes 
referred to as eccentric anomaly and is related to the polar angle 8 
by the following equation: 

btan<j> = atan(#— a). (25) 



3.3.3 Fitting the limb and the terminator 

Note that the the limb and terminator can be represented in a single 
parametric form [iirtf! 1 ), 2/lt(<^)] = /lt(^) using the equations 
given in section|2] This shape can be fitted by considering the three 
following transformations: homothetic transformation with scale 
factor c, rotation with angle a and translation by (x c , Vc) 



f(<f>; [x c ,Va,c,a]) = 
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In this case the global parameters for f LT (4>) are the geometric 
parameters /3®, /3q, AA, r e , r v introduced in section[2](as well as 
the distance from the observer to the planet to convert projected 
length into pixels units). These parameters may be determined by 
e.g. SPICE. The unknown global parameters to optimise for f(4>) 
are the planet centre (x c ,y c ), the scale factor c and tilt angle a. 



3.4 Algorithm 

Phase (i) is performed once while phases (ii) and (iii) can be iter- 
ated. The iteration process improves the selection of seeds for the fit 
and removes any "outliers" that might be included using the crude 
estimates for the parameters, therefore improving the accuracy of 
the fit. The tolerance on fractional improvement of Q defined in 
Eq. \22\ is set to 10 -3 . Such tolerance ensures convergence of the 
Levenberg-Marquardt algorithm in a few iterations leading to an 
accuracy of the estimate of the centre coordinates and the radii of 
the order of one pixel or better. 



4 APPLICATION TO PLANETARY IMAGES 
4.1 Ellipse for complete planet 

The image presented here (see upper panel of Fig. |2| to illustrate 
the location method of a complete planetary disc has been obtained 
using the 3.8 m UKIRT at Mauna Kea observatory, Hawaii, with 
the near-IR UIST guide camera (Ra msay Howat et al.|2004j l. This 
image has not been flux calibrated but the sky background noise 
has been subtracted, and the intensities are thus in arbitrary units. 

UIST is a 1-5 /J,m imager-spectrometer with a 
1024x1024 pixels InSb array. In imaging mode there are 
two plate scales available, with resolution 0.12 arc sec/pixel and 
0.06 arc sec pixel -1 , giving fields of view of 2x2arcmin 2 and 
lxl arc min 2 respectively. 

UIST was used to observe Jupiter at a resolution of 
0.12 arc sec pixel -1 , with the Brackett alpha filter (50 percent 
cut-on at 4.024 /on and 50 percent cut-off at 4.078 (jm) in expo- 
sures of 10 s. The Brackett line is an IR emission line of the H atom. 
Thus this emission should contribute many of the photons as well as 
H J . In this part of the IR spectrum, the emission of the giant planets 
is dominated by several lines of , and the spectral measurement 
of individual lines allows determination of temperatures and 
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Original image 




-300 -200 -100 100 200 300 

Figure 2. In the upper panel, median-filtered IR image collected by UKIRT 
(in arbitrary intensity unit) and in the lower panel the resulting Voronoi 
tessellation from the VOISE algorithm. The 1612 tiles of the tessellation 
are uniformly coloured using the median intensity of the pixels that are 
lying within each polygon. The axes are labelled in pixels unit and the point 
with coordinates (0, 0) is the centre of the image. The colour code for both 
images is shown in Fig.R] 



column densities of the planet (Miller et al. 20061. The UIST cam- 
era was used in conjunction with the dual-beam polarimeter module 
IRPOL2 for spectropolarimetry measurements under an observa- 
tion campaign of Jupiter on August 4, 2008. 

Fig. [2] illustrates phase (i) of detection of the limb using 
VOISE on an image collected by UKIRT at 10:13:00 UT. The size 
of the image is 679x639 pixels and it has been pre-processed by 
a nonlinear filter — a median filter — of size 11x11 pixels in order 
to lower noise in the image ( Gonzalez & Woods 2007 1. Whenever 
such noise filter is used as pre-processing to the VOISE segmen- 



L =16 C(0.0,0.0) a=289.9 b=271.1 a=0 £(0.90,1.10) l 




Figure 3. Selection of the seeds from the VOISE tessellation for two it- 
erations of the phases (ii) and (iii). During the first iteration (upper plot) 
666 seeds are selected as neighbours of the limb while for the second iter- 
ation the numbers of seeds considered for the fit is reduced to 583 seeds. 
The size of each coloured marker is proportional to the surface area of the 
selected polygon. The limits of the torus are shown in thick red lines while 
the nominal ellipse is shown as a thin red line and the red cross is the centre 
of the torus. 



tation, the size of the mask should be chosen to be larger than the 
minimum seed distance d m to be of any effect. The main idea of 
this filter is to slide a window with specified size and replace each 
centre pixel of the window by the median of the pixels lying in 
the w indow. The VOISE parameters iGuio. R and Achilleos, N. 
2009j > are (i) division phase: d 2 m = 9 pixels , pd = 97 percent 
(ii) merging phase: pm = 50 percent, A/j, = 20 percent and 
AH = 30 percent (iii) two iterations in the regularisation phase. 
The resulting segmentation contains 1612 polygons (lower panel 
in Fig. [2}. Note the compactness of the polygons in regions with 
small length scales along the limb and near the equator. 

Fig. [3] illustrates phase (ii) of selecting the set of points from 
the segmentation to be used as the neighbourhood of the plane- 
tary limb. The upper panel shows the first iteration, using crude 
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Table 1. Parameters of the fitted ellipse resulting from two iterations of phases (ii) and (iii). The ellipses are shown in Fig. [4] Note that the tilt angle has not 
been fitted and has been fixed to a = 0. # iter is the number of iterations performed in order to converge with the prescribed tolerance and the normalised x 2 
provides an indication of the goodness of the fit. The "guess", "fit 1" and "fit 2" ellipses are shown in Fig. [4] 
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Figure 4. Initial guess for the ellipse and the ellipses resulting from the two 
successive fits for the original image shown in Fig. [2] The points selected 
for each fit are as shown in Fig. [5] 
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Figure 5. Image scaled in arc sec and the projection in the sky-plane of a 
latitude-longitude grid of Jupiter computed using data from SPICE and the 
fitted parameters of Jupiter's disc. 



estimates of the ellipse parameters to define the large torus with 
£m = 0.9 and em = 1.1, and a relatively large scale length param- 
eter Cm = 16 pixels. The lower panel shows the selection process 
for the second iteration with parameters provided by the result of 
the first fit. The torus has been re-centred and its thickness reduced 
by setting e m = 0.9 and em = 1.05. The scale length parame- 
ter has also been reduced to Cm = 12 pixels (to be compared to 
d m = 3 pixels). In both iterations, seeds in the neighbourhood of 
the faint emission outside the limb near the South pole have been 
rejected whenever the polar angle of the seed ip with respect to the 
centre C(x c ,y c ) is in the range — 80°<</?< — 65°, i.e. the seed 
lies inside the grey shaded sector depicted in Fig. [3] 

Fig. [4] illustrates phase (iii) consisting of the nonlinear fitting 
of the selected points (shown in Fig. [3} to an ellipse in parametric 
form given by Eq. l |24| (. The planetary disc modelled as a single el- 
lipse is justified in the situation where the disc is nearly fully illumi- 
nated, as is the case here. The curve labelled "guess" corresponds to 
crude estimates of the ellipse parameters, i.e. the coordinates of the 
planet centre correspond to the centre of the image and the equato- 
rial and polar radii are derived using SPICE. The curve "fitl" is the 
curve with parameters after first fit, i.e. with the seeds as seen in the 
upper panel of Fig.[3]and the the curve labelled "fit2" corresponds 
to seeds as seen in the lower panel of Fig. [3] 

The parameters, error estimates and fitting parameters are 
given in Table[T] It is interesting to note that the estimated param- 
eters related to the x-direction (x c and a) have smaller errors com- 
pared to the parameters relates to the y-direction (y c and b) which 
is a consequence of the large sampling of seeds around the equator. 
We have also checked, for consistency, that the curve parameters 
{4>i\ and the global parameters p together with their errors {A<fii} 



and Ap provide error in positioning of the points consistent with 
the values used for the weights rjj in Eq. j22}. 

Fig. [5] presents the image with a (planetocentric) latitude- 
longitude grid (with 10 ° step in latitude and 20 ° in longitude) 
computed using the coordinates of the centre of Jupiter's disc ob- 
tained from the nonlinear fitting and the projection geometry com- 
puted using SPICE. The Central Meridian Longitude (CML) of 
Jupiter at the time of the observation is 157.5 deg. The limb is 
shown on the right side of the planet (red solid line) while the ter- 
minator is shown on the left side (green solid line). The parameters 
are /3© = -1.5°, fig = -1.4° and AA = 5.3°. The equato- 
rial radius for the final fit is r e = 70958 km with an eccentricity 
e = 0.32 while the values provided by SPICE are r e — 71492 km 
and e = 0.35. The illuminated limb in the considered near-IR 
waveband is thus slightly smaller than Jupiter's 1 Bar pressure sur- 
face as given by SPICE. 

4.2 Circle for partial planet 

The image presented in this section (upper panel in Fig. [6} has been 
chosen to illustrate the case with partial occlusion of the planetary 
disc. It was collected with NASA's 3.8 m IRTF at Mauna Kea ob- 
servatory, Hawaii using the imaging facility ( |Shure et al.| 19 941 at 
wavelength 3.43 /im (a wavelength sensitive to H3+). The image 
was collected during a campaign on June 28, 1995 at 11:14:52 UT. 
This image has not been flux calibrated but the sky background 
noise has been subtracted, and the intensities are thus in arbitrary 
units. 

The NSFCAM is a 1-5 imager with a 256x256 pixels 
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Original image 




Figure 6. In the top panel, IRTF image processed by a median filter fol- 
lowed by histogram equalisation (in arbitrary intensity unit). In the lower 
panel the result of the segmentation by VOISE. The colour code for the 
494 tiles is the same as for the segmentation presented in the lower panel 
in Fig. [2] The axes are labelled in pixels unit and coordinates (0, 0) cor- 
respond to the centre of the planet provided by the original Image Reduc- 
tion and Analysis Facility (IRAF) data reduction, taking into account an 
estimate based on the telescope pointing (Satoh & Connerney 19991. The 
colour code for both images is shown in Fig. [8] 



InSb detector. Three different magnifications are avail- 
able: 0.3arcsec pixel -1 , 0.15 arc sec pixel -1 and 
0.06 arc sec pixel -1 corresponding to a field of view of 
76.8 arc sec, 37.9 arc sec and 14.1 arc sec respectively. 
The NSFCam has been upgraded (NSFCam 2) with a 



Table 2. Resulting parameters for the circle during the two iterations of fit. 
The number of iterations and the resulting normalised \' 2 are also given. 
The corresponding "guess", "fit 1" and "fit 2" circles are shown in Fig. [8] 







Vc 


R 


# iter 


x 2 


guess 


0.0 


0.0 


143.9 






fit 1 


-0.3±0.9 


2.5±0.9 


149.8±0.9 


3 


2.87 


fit 2 


-0.3±0.8 


3.5±0.8 


150.3±0.8 


3 


1.71 



2048 x 2048 pixels Hawaii-2RG detector. The image scale 
will be 0.04 arc sec pixel -1 with field of view 80x80 arc sec 2 . 

Fig. [6] shows the results of phase (i) of the method using 
VOISE on an image collected by UKIRT at 0721 UT. The im- 
age has size 256 x 256 pixels. Note that the image has been pre- 
processed by a median filter of size 7x7 pixels to lower noise level, 
followed by a histogram equalisation (Gonzalez & Woods 2007| >. 
The histogram equalisation is performed in order to increase the 
global contrast of the original image (which is shown in Fig. [9}. 
It consists of a nonlinear adjustment of the intensities in order to 
better distribute the image intensity histogram and accomplishes 
this by effectively "spreading out" the most frequent intensity val- 
ues. Alternatively, the contrast in the low intensity range can also 
be enhanced by taking the logarithm of the ratio of the image pix- 
els relative to the estimated noise level, if available. Note that we 
haven't pre-processed the UKIRT image for the first example as 
the limb boundary was already substantially more intense than the 
background. The VOISE parameters have been set to (i) division 
phase: d 2 n — 4 pixels 2 , pr> = 98 percent (ii) merging phase: 
Pm = 50 percent, A/i = 20 percent and AH = 30 percent 
(iii) two iterations in the regularisation phase. 

Fig. [7] illustrates phase (ii) of selecting the set of points from 
the Voronoi map to be used as within the limb neighbourhood. Note 
that the points with polar angle <p with respect to the centre estimate 
C(x c ,y c ) such that — 115 <<,5< — 65° (inside the grey shaded 
sector in Fig. [7]( are filtered out to avoid bias from the seeds cor- 
responding to the the emission outside the limb which has been 
highlighted by the histogram equalisation. 

Fig. [8] illustrates phase (iii) consisting of the nonlinear fitting 
of the selected points (shown in Fig. ^ to a circle in parametric 
form Eq. J23[ >. The global parameters of the circle fitted together 
with estimates for the error are given in Table [2] The circle as a 
model for the disc is justified in situations where only a portion 
of the planetary disc is in the field of view, as it is the case in the 
present image. Note also that in this case the distribution of the 
seeds is uniform from the equator to the South pole and therefore 
the estimated parameters related to the a;-direction have similar er- 
rors as those related to the {/-direction. 

Fig. [9] presents the image with the same (planetocentric) 
latitude-longitude grid resolution as in Fig. [5] calculated with the 
coordinates of the centre and radius of Jupiter's disc obtained from 
the nonlinear fitting, and the projection geometry from SPICE. The 
CML of Jupiter at the time of the observation is 13.5 deg. The limb 
and terminator are shown as red and green solid lines respectively, 
and the thick black line is the noon-meridian. The parameters pro- 
vided by SPICE are /3 ffi = -2.9 °, /3 = -2.8 and AA = 5.3 °. 
The radius for the final fit ("fit 2") is r e — 72466 km. The illumi- 
nated limb in this waveband is slightly larger than Jupiter's 1 bar 
pressure surface. We also tried to fit an ellipse and it leads to a simi- 
lar estimate for the centre coordinates but with larger error bars due 
to the increased degrees of freedom. 
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L =10 C(0.0,0.0) R=143.9 £(0.90,1.10) 
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Figure 7. Seeds selection for two iterations of phase (ii). For details see 
Fig.0 

5 DISCUSSION 

We have presented a novel semi-automatic method to estimate ac- 
curately, and objectively, the disc parameters in an image of an 
illuminated planetary disc. The method is based on the "best" fit 
of a set of points selected from a segmentation map generated by 
VOISE to a curve described in a parametric form. 

The segmentation phase can be improved by pre-processing 
the image using different techniques such as noise filtering and con- 
trast adjustment. 

Basic shapes to describe the boundary of a planetary disc in- 
clude the circle and the ellipse. We also provide analytic expres- 
sions for the projection in the sky-plane of the limb and terminator 
of a planet modelled as ellipsoid. These expressions can easily be 
used to describe both limb and terminator as single curve in para- 
metric form. 

Note that the VOISE algorithm generates "intermediate" tes- 
sellations, one at the end of the division phase and one at the end 




-50 50 100 150 
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Figure 8. Initial guess for the circle and the circles "fit 1" and "fit 2" result- 
ing from two iterations of the phases (ii) and (iii). The selected points for 
each iteration are shown in Fig. [7] 
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Figure 9. Original image scaled in arc sec together with the sky-plane pro- 
jection of a latitude-longitude grid computed using data from SPICE and 
the fitted parameters of Jupiter's disc, i.e. the circle parameters. 



of the merging phase. It is worth noting that fitting an ellipse gives 
the best result (smallest x 2 ) f° r the regularised tessellation (i.e. af- 
ter merging), but errors in the fitted parameters are smaller when 
considering the map at the end of the division phase. The largest 
X 2 is obtained for the map generated at the end of the merging 
phase. This confirms that the tessellation obtained after division is 
optimum for our purposes. The reason for this is that VOISE merg- 
ing generates more regular polygons, but very slightly degrades the 
position information from the division phase. 

We have shown that our novel objective method to locate the 
planetary disc on images provides improved estimates of the centre 
position (as compared to the guide star catalogue) as well as the 
altitude when the disc is illuminated for the corresponding obser- 
vational waveband. 
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We also showed that the use of histogram equalisation en- 
hances the auroral emission outside the limb and therefore allows 
a better and unbiased estimate of the limb by allowing removal of 
points from this auroral emission region. 

The software implementing this method is written in MAT- 
LAB® and can be made available by request to the authors. 
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